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Abstract 

We review a molecular dynamics method for nucleon many-body systems 
called the quantum molecular dynamics (QMD) and our studies using this method. 
These studies address the structure and the dynamics of nuclear matter relevant 
to the neutron star crusts, supernova cores, and heavy-ion collisions. A key ad- 
vantage of QMD is that we can study dynamical processes of nucleon many-body 
systems without any assumptions on the nuclear structure. First we focus on 
the inhomogeneous structures of low-density nuclear matter consisting not only 
of spherical nuclei but also of nuclear "pasta", i.e., rod-like and slab-like nuclei. 
We show that the pasta phases can appear in the ground and equilibrium states 
of nuclear matter without assuming nuclear shape. Next we show our simulation 
of compression of nuclear matter which corresponds to the collapsing stage of su- 
pernovae. With increase of density, a crystalline solid of spherical nuclei change 
to a triangular lattice of rods by connecting neighboring nuclei. Finally, we dis- 
cuss the fragment formation in expanding nuclear matter. Our results suggest 
that a generally accepted scenario based on the liquid-gas phase transition is not 
plausible at lower temperatures. 

§1. Introduction 

Due to the progress of computers, numerical simulations became increasingly 
capable in tackling complicated problems in nuclear physics. Generally, numerical 
simulations can be classified into two types: macroscopic and microscopic simu- 
lations. The former, macroscopic simulations, deal directly with the macroscopic 
quantities which we are interested in. We need to introduce physics models which 
describe how the quantities are connected with each other. On the other hand, 
microscopic simulations are based on the degrees of freedom of the constituent el- 
ements. The necessary inputs are equations of motion and the interactions among 
the elements. The properties of the total system are obtained later by analyzing 
the resultant information of these constituent elements. Microscopic simulations 
have several advantages: 1) We need only a few assumptions on the model. 2) 
We may obtain unexpected results. 3) We may find a physical principle (the law 
governing the elements) if we obtain suitable observables. 
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This review article is about the molecular dynamics (MD) simulations of nu- 
clear matter. First, we give an overview of the history of the MD models in nuclear 
physics. The simulation study of nuclear dynamics originated from the formula- 
tion of the time-dependent Hartree-Fock (TDHF) theory in 1930.^ TDHF deals 
with the time evolution of many-fermion systems and is an approximation of the 
time-dependent Schrodinger equation with the use of a single Slater determinant. 
However, it was in the 1970's that the TDHF was first solved numerically.'^l' Due 
to the limitation of computer power, only low-energy phenomena were studied 
in the early stage. As the computational power drastically increased after the 
1980's, applications to higher-energy phenomena with larger numbers of degrees 
of freedom and also improvements of the framework to include correlations, etc., 
have been made. However, TDHF cannot describe the heavy-ion collision process 
at higher energies where the degrees of freedom drastically increase. The reaction 
mechanism depends on the incident energy and the impact parameter, as shown 
in Fig. [TJ Particularly, the Fermi energy is the key quantity to characterize the 
mechanism. It is because the reaction mechanism is determined by the competi- 
tion between Fermi motion inside the nuclei and the relative motion of colliding 
nuclei. Above the Fermi energy (medium - high energy), the region where the 
colliding nuclei overlap with each other will break into fragments, i.e., fragmen- 
tation occurs. This fragmentation process is driven by the energy of the relative 
motion between the two nuclei converted into thermal excitations, which are gen- 
erated by two-body collisions. Each collision is regarded as a transition from a 
Slater determinant into another. Such a change of the wavefunction cannot be 
described by TDHF with a single Slater determinant. 
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Fig. 1. Schematic diagram of heavy-ion reaction mechanism. 

At higher energies, the above-mentioned two-particle collision becomes an im- 
portant process which determines the reaction mechanism in the heavy-ion col- 
lisions. The time evolution of the phase-space distribution function in heavy-ion 
collisions is described by a Boltzmann-type equation of motion (EOM), with a 
smooth change by a Newtonian equation and dissipation by the two-body collision 
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process, which is called Boltzmann-Uehling-Uhlenbeck (BUU), Vlasov-Uehling- 
Uhlenbeck (VUU) , or Boltzmann-Nordheim-Vlasov (BNV) equation. If one omits 
the collision term, this framework can be regarded as a classical limit (ft — > 0) of 
TDHF,'2J i.e., the Vlasov equation. This Boltzmann-type equation can be numeri- 
cally solved by a test-particle method: The fluid elements in a 6-dimensional phase 
space are replaced by a classical particle and the phase-space distribution function 
is obtained by counting the number of those particles in the 6-dimensional mesh 
(three dimensions for coordinate space and the other three dimensions for mo- 
mentum space). For sufficiently large numbers of test particles, the 6-dimensional 
particle density is conserved throughout the time evolution in a mean-field po- 
tential described by the Vlasov equation. The two-body collision process of test 
particles violates the conservation of 6-dimensional phase-space distribution. In 
the 1980's many works on the heavy-ion collisions in the medium - high energy 
region have been made via BUU simulations with the test-particle method. 




Fig. 2. Schematic explanation on the time evolution of two-particle correlation. The curves indicate 
the distribution functions and the circles are their representative test particles. Refer to the 
text for details. 



Molecular dynamics simulation for nuclear systems has been developed in the 
late 1980's. Aichelin and Stocker have proposed quantum molecular dynamics 
(QMD) model to simulate heavy-ion collisions from medium to high energies.'^ 
This framework is obtained by reducing the number of test particles in BUU 
simulations so that each particle represents one nucleon. By this reduction of 
test particles, it became possible to describe many-body correlation of the system. 
Let us take the two-body correlation as an example. If two nucleons stay within 
a distance so that their distribution functions (solid curves in Fig. [2]) overlap 
with each other, the representative two test particles (circles in Fig. [2]) can be 
very close to each other (corresponding to the left panel of Fig. ^ and also can 
be far (right panel). In the former case, those two nucleons may be bound to 
form a cluster while, in the latter case, they may diverge from each other as 
time passes. However, with a huge number of test particles, we obtain only one 
time evolution of the distribution function which corresponds to the average of 
many events. This example of the two-body correlation effect on the fragment 
formation is related to the variation between time evolutions of different events. 
It is also natural that the two-body correlation is important for the description of 
spatial fluctuation in a event. In the mean-field calculation, on the other hand, 
both fluctuations in space and the variation between events will be washed out. 

The QMD model assumes a direct product of nucleon single-particle wave- 



functions as a total wavefunction and a Lagrangian with a non-relativistic kinetic 
energy and a potential energy from effective interactions among nucleons. The 
single-particle wavefunction is assumed to be a Gaussian wavepacket with a fixed 
width. The EOM of the wavefunction is derived from a variational principle with 
the above Lagrangian, and results in a classical EOM with a Hamiltonian as a 
function of coordinates of those Gaussian wavepackets. In the QMD model, the 
stochastic two-body collision process is added to the time evolution by the Hamil- 
ton EOM. The final state of the two-body collision process is checked so that it 
obeys the Pauli principle, i.e., the condition on the phase-space density. 

QMD is named "quantum" due to (1) the many-body correlation or fluctua- 
tion in density caused by the EOM and the collision term, (2) the stochasticity 
in the collision process, (3) the Pauli blocking in the final state of collision, and 
(4) the use of Gaussian wavepackets for single-particle wavefunctions. However, 
the actual feature of QMD simulation is rather classical. First, the time evolu- 
tion of the system concerns only the centroids of wavepackets. Their width in 
the coordinate space, which is a fixed parameter, appears only in the interaction 
among the particles by means of the double folding. The width in the momentum 
space gives rise to a part of the kinetic energy. However, this energy is spurious, 
i.e., it will never be effective, since it is constant during the time evolution. Sec- 
ond, the Pauli principle, which yields the fermionic momentum distribution, is 
not basically taken into account. When the number of the test particles per nu- 
cleon is large enough in the Boltzmann type simulation, the phase-space density 
is conserved in the moving frame of a fluid element. In spite of the many-body 
nature obtained by the reduction of the number of test particles, it sacrifices the 
fermionic nature of the system. 

One of the most serious problems of QMD is in the description of the ground 
state. Due to the lack of fermionic characteristics, the energy minimum states 
of QMD model violate the Pauli principle, and all the particles degenerate into 
zero in the momentum space so that they overestimate the binding energy. We 
cannot use the energy-minimum states as the initial conditions of collision sim- 
ulations. If we prepare initial conditions with appropriate binding energies, the 
constituent nucleons would be moving. This motion makes the initial condition 
unstable against the emission of nucleons. To take the fermionic characteris- 
tics into account, we need to introduce explicitly the antisymmetrization of the 
wavefunction.EJ'lSJS 

Fermionic molecular dynamics (FMD)^' and antisymmetrized molecular dy- 
namics (AMDJS) have been proposed in 1990 and 1992, respectively. They assume 
a Slater determinant of Gaussian wavepackets as the wavefunction of the system. 
In FMD, the widths of nucleons are time-dependent variables and the kinetic en- 
ergies of wavepackets are not spurious. In AMD, the widths of wavepackets are 
constant in time but the zero-point center-of-mass kinetic energies of fragments 
are removed in a phenomenological way. They have succeeded in describing the 
ground state properties of light nuclei as well as the dynamical processes of low- 
energy heavy-ion collisions. The problem is, however, a huge amount of computing 
cost to solve the equations of motion of FMD and AMD, which is proportional 
to the fourth power of the particle number N (cf. oc N'^ for QMD). Thus the use 
of FMD and AMD has been limited to small systems with the total number of 
particles up to a few hundreds. 

In this situation, a new phenomenological way to mimic the Pauli principle 
was introduced in QMD.I^ Wilets et al^ and then Dorso et al^^ developed 



a repulsive two-body potential so-called the Pauli potential. It is a function of 
not only the distance in the coordinate space, but also of the distance in the 
momentum space. This repulsive potential acts between nucleons with the same 
spin and isospin so that it prevents those particles from coming close in the phase 
space. Note that, in this framework, simulated ideal Fermi gases contain the 
potential energy which comes from the Pauli potential. It is counted as a part 
of the nuclear potential energy when one determines the parameters of effective 
potential. Due to the momentum dependence of the Pauli potential, constituent 
nucleons have non-zero values of the momentum in the ground state keeping 
their velocities at zero; thus the above-mentioned spurious emission of nucleons 
is avoided. 

Since the appearance of the QMD model with the Pauli potential, it became 
possible to carry out simulations of systems with a large number of nucleons. One 
interesting target was low-density (below the saturation density) nuclear matter 
in compact stars such as crusts of neutron stars and cores of supernovae. In 
low-density nuclear matter, exotic structures called "nuclear pasta" have been 
predicted by Ravenhall et al^^ and Hashimoto et al^^ There, nuclear matter 
cannot be uniform due to a negative partial pressure of nucleons and should 
be clusterized. With increasing density, the shape of the cluster changes from 
droplet, rod, slab, tube, bubble, and then uniform. The name of "pasta" comes 
from the similarity of the rods to "spaghetti" and the slabs to "lasagna" , etc. Since 
Ravenhall et al. and Hashimoto et al. have proposed, many works have been done 
on nuclear matter with pasta structures. Most of them are based on the Wigner- 
Seitz (WS) approximation, in which a unit cell with the dimensionality 1,2, and 
3 is replaced by the same volume of the plate, cylinder, and sphere, respectively. 
The WS approximation is useful and saves much CPU time. However, the use of 
WS cell should be a strong constraint on the structure and only simple structures 
are allowed. On the other hand, MD simulation is a microscopic framework which 
does not need any assumption on the structure and the reaction mechanism. Since 
a QMD model is capable of simulating systems with a huge number of particles, we 
have applied it to low-density nuclear matter and its inhomogeneous structures. 
In Sec. [3l we present the results of our QMD simulations. 

Apart from adapting QMD to low-energy phenomena, other efforts have been 
made to the opposite direction, i.e., an attempt to describe high-energy phe- 
nomena. Sorge et al. have proposed relativistic quantum molecular dynamics 
(RQMD^Sli in 1989. The main improvements of RQMD from QMD are (1) the 
Lorentz covariance in the interaction, kinematics, and the two-body collisions and 
(2) the inclusion of baryon resonances, strange particles, and the string excitations 
in the two-body collision process. Simulations of high-energy heavy-ion collisions 
have been carried out to analyze experiments with E/Ak,\- 200 GeV at the SIS, 
AGS, and SPS facilities. At further high energies, interaction between particles 
becomes less important and the production of mesons and excitations of baryon 
resonances are essentially important. A set of computational codes called "ultra- 
relativistic QMD" (UrQMD) is developed with a collision term highly tuned-up 
to include various kinds of baryons, mesons, and their excited states.l^ It is 
distributed on the internet and is often used by many people for simulations of 
heavy-ion collisions at RHIC experiments. 



§2. Molecular dynamics approach to nuclear matter 



2.1. The total wavefunction and the equation of motion 

In QMD, each nucleon state is represented by a Gaussian wavefunction of 
width A, 
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where R; and P; are the centers of position and momentum of i th nucleon, respec- 
tively. The total wavefunction is assumed to be a direct product of these wave- 
functions. Thus the one-body distribution function is obtained by the Wigner 
transform of the wavefunction, 
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The equations of motion of R^ and Pi are given by the Hamiltonian equations 
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and a stochastic nucleon-nucleon collision term. Hamiltonian H consists of the 
kinetic energy and the energy of the two-body effective interactions. 

2.2. Effective interactions 

Our Hamiltonian is separated into several parts as follows, 

H = T+ V^Pauli + V^ocal + VmD , (5) 

where T, VpauU, Mocai, and Vmd are the kinetic energy, the Pauli potential, the lo- 
cal (momentum-independent) potential, and the momentum-dependent potential 
parts, respectively. 

The Pauli potentiaP^''i^''i^''i^ is introduced to mimic the fermionic proper- 
ties in a semiclassical way. This phenomenological potential prohibits nucleons 
of the same spin a and isospin r from coming close to each other in the phase 
space. Here we employ the Gaussian form of the Pauli potentiaP^'^i^ as 



Vpauli = TT^P 
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In the local potential part, we adopt the Skyrme type with the Coulomb and 
the symmetry terms as explained in Eq. (5) of Ref. [T^ . 
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In the above equation, po is the normal nuclear density (~ 0.165fm^'^), Ci is 1 for 
protons and for neutrons, while (pi) and (pi) are overlaps of density with other 
nucleons defined as 
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It is known that the nucleon-nucleon interaction has a strong momentum- 
dependence [see Fig. [3]- We have chosen the form of the momentum-dependent 
term as a Fock term of the Yukawa-type interaction. We divide this interaction 
into two ranges so as to fit the eff'ective mass and the energy dependence of the 
real part of the optical potential, as 
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Parametrization of the constants in the above effective interactions will be 
discussed in Sec. 12.51 

2.3. Energy minimum state 

For the Hamiltonian with the Pauli potential, we can define the ground state 
as an energy-minimum state of the system. To get the energy-minimum configu- 
ration, we use the following damping equations of motion, 

^ _ dH dH ■ _ dH dH 

R.-^-MR^, P.--^-MP^, (12) 

where /iR and /ip arc the damping coefficients with positive values. 

We first distribute the particles randomly in the phase space and cool down 
the system according to the damping equations of motion until the energy reaches 
the minimum value. Sometimes the system is trapped in a local minimum. We 
thus try again and again this cooling procedure with a different initial state and 
seek the global energy minimum state. 

For finite nuclei and nuclear matter above the saturation density, this proce- 
dure works well. For nuclear matter at subsaturation densities, however, there 
are many local minimum states around the true ground state, which differ from 
the ground state in the details of the surface configuration of clusters. Since the 
energy difference from the ground state is the order of 10 keV/nucleon in this 
case, we accept these states as ground states and neglect the small differences of 
the configuration. 



2.4. Periodic boundary conditions 

In order to simulate infinite nuclear matter with finite numbers of particles, 
we use a cubic cell with periodic boundary conditions. The size of the cell is 
determined from the average density and the particle number. The periodic 
boundary conditions can be introduced as follows: We prepare 26 (= 3'^ — 1) 
surrounding cells, which are copies of the central cell. The particles in the central 
cell move according to the interaction with all particles in the same cell and in 
the surrounding cells as well. The particles in the surrounding cells obey exactly 
the same motions as those in the central cell. Thus the Hamiltonian per cell is 
written as 



H = 
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where Ti is one-body part (kinetic energy) , H^j ' is the two body part of the 
Hamiltonian and Lcoii are the relative position of surrounding cells from the cen- 
ter. Note that the indices "cell" runs from (the central cell) to 26 (surrounding 
cells) and Lq = 0. 

2.5. Parametrization of the constants 

We have twelve parameters in the effective interactions of the Hamiltonian 
([5]), i.e., Cp,(7oiPOi CK,P,T,Cs, Cii'', Ccx\ A^i, /^2, and the Gaussian width A. We 
fix these constants to reproduce properties of the ground state of finite nuclei and 
saturation properties of nuclear matter. 

We first determine the parameters, qo-,po-, and Cp, of Pauli potential apart 
from the other effective interactions, by fitting the kinetic energy of simulated 
matter to the energy of the ideal Fermi gas at zero temperature and at various 
densities. For this, we define the free Fermi gas system as a ground state for 
the Hamiltonian including only the kinetic energy and the Pauli potential by 
making use of the damping equations of motion ([T2|) and the periodic boundary 
conditions with 1024 particles in a cell. In Fig. |4j we show the kinetic energy 
(the solid squares) and the total energy (the open squares) obtained by the Pauli 
potential with 



Cp = 207 MeV, po = 120 MeV, and qa = 1.644 fm. 
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In the same figure, we plot the energy of the Fermi gas with a solid line. Although 
there are some other parameter sets which can reproduce the ideal energies of the 
Fermi gas using the same form of the Pauli potential, i.e., that used in Ref. |5]), 
we choose the above parameter set to get good properties of the ground state 
of finite nuclei with other effective interaction terms particularly in combination 
with the momentum-dependent interaction. 

Among remaining nine conditions, four are for the momentum-dependent in- 
teraction as follows. We calculate the single particle potential of momentum p in 
nuclear matter at the normal nuclear density, which leads to 
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Fig. 3. Momentum dependence of the potential energy of 
experimental data and the present QMD model. This 
figure is taken from Ref. I20p . 



P/Po 

4. Density dependence of the en- 
ergy per particle of an ideal Fermi 
gas. The solid curve shows the ex- 
act value and the symbols show the 
QMD results using the Hamiltonian 
which contains only the kinetic en- 
ergy and the Pauli potential. The 
filled squares show the kinetic en- 
ergy and the open squares show the 
total energy. This figure is taken 
from Ref.l^. 
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We fit the energy dependence of this potential to experimental data. In Fig. [31 
we plot the energy dependence of the real part of the optical potential (the open 
circles and squares) obtained from the experimental data of Hama et al^^ for 
proton-nucleus elastic scattering. To fit the data, we impose three constraints, 
i.e., t/(0) = -80 MeV, U{p) = at Eux> = 200 MeV, and U{p oo) = a + jS = 77 
MeV. For another condition, we take the effective mass m* = 0.8 m at p — po, 
where 

1 1 /lac/MD 



p dp 



(17) 



Other three conditions are coming from the saturation condition, i.e., the 
energy per nucleon E/A = —16 MeV and the incompressibility K at p ^ po. 

There are two parameters left. One is the symmetry energy coefficient Cg, 
which we take 25 MeV to get a reasonable value of the symmetry energy 34.6 
MeV for nuclear matter at the saturation density. The other is the width of the 
Gaussian wavepacket A, which is chosen to get appropriate ground state properties 
of finite nuclei and infinite nuclear matter below the saturation density without 
changing those of uniform nuclear matter above the saturation density. We then 
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choose these parameters to give a good fitting to the binding energies of finite 
nuclei plotted in Fig. [S] 




Fig. 5. Binding energies of nuclei calculated with our QMD model. There are three sets of param- 
eters, soft, medium, and hard, which refer to the stiffness of nuclear matter at the saturation. 
This figure is taken from Ref. I20|l. 
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Fig. 6. Energy per nucleon calculated with our QMD model. This figure is taken from Ref. I20|l. 



It should be noted here that we cannot determine these parameters from the 
above conditions in an analytical way, since the Fermi distribution is not exactly 
achieved by the Pauli potential and the additional potential energy included in 
the Pauli potential. Thus we simulate nuclear matter by QMD with the periodic 
boundary conditions using 1024 particles in a cell. We search the energy minimum 
state by the damping equations of motion ([T^ as discussed above and adjust the 
parameters. By this method, we have fixed three parameter sets corresponding 
to equations of state (EOSs) with different values of incompressibility K: Soft 



(K=210 MeV), Medium (if =280 MeV), and Hard {K=380 MeV) EOSs.iSSJ These 
values of the incompressibihty K are extracted from the curvature of the energy 
per nucleon (shown in Fig. ^ at the saturation density by fitting to the following 
parabolic form, 

E/A - (P - Pof - 16 MeV . (18) 

The single particle potential shown in Fig. |3]are also calculated by the simu- 
lated nuclear matter with the Pauli potential with the effective interactions ([5]). 
The results are denoted by the crosses in Fig. [3] and agree well with the single 
particle potential given by Eq. (jlSp except for the low energy part, where the 
Pauli potential is effective. Though this result in Fig. [3] is obtained with the 
parameter set of Medium EOS, results with Soft and Hard EOSs are the same as 
Medium EOS within 2 MeV for the whole energy region. 

In Fig. [SI we plot the binding energies of the ground state of finite nuclei 
obtained by the damping equations of motion (jl2p with three parameter sets, 
i.e., Soft (the long dashed line). Medium (the dashed line), and Hard (the solid 
line) EOSs. All of them describe well the global trend of the binding energies of 
various nuclei except for light nuclei from ^^C to ^"Ne. The disagreement for the 
light nuclei might be due to the individual structures (shell structures, cluster 
structures, etc.) of these light nuclei, which are not well described by the present 
QMD. The parameters for Medium EOS which we use for our QMD calculations 
are listed in Table HI 



Table I. Effective interaction parameter set (7^=280 MeV). 



a (MeV) 


-92.86 


P (MeV) 


169.28 


7 


4/3 


Ceo (MeV) 


25.0 


(MeV) 


-258.5 


C,^(2) (MeV) 


375.6 


pi (MeV) 


2.35 


P2 (MeV) 


0.4 


Cp (MeV) 


115.0 






po (MeV) 


207.0 


qo (fm) 


1.644 


A (fm) 


1.45 







2.5.1. Electron background 

The total charge of the system should be zero, i.e., the numbers of electrons 
and protons in the cell are equal. As the distance between nuclei is small compared 
with the Thomas- Fermi screening length by electrons A^p, where 



4^A^/ = f^, (19) 
op, 



the distribution of electrons should be almost uniform. In fact, the inhomogeneity 
of electron distribution is found to be small by our studied^ which include 
the screening by electrons. 

Therefore, it is natural to treat electrons as a uniform background. In this case, 
we should take account of the long-range contributions of the Coulomb interaction 
between protons. However, for simplicity a cutoff of the Coulomb interaction, e.g., 
in the form of screened Coulomb interaction, is sometimes used .'^''^''23 In Sec. 
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13.11 we employ this "screened" Coulomb interaction among protons, 

^c" = Y E // rf^^.^V /""^',';"':'^''^'^"^ P.(r.)ft(r.) , (20) 

where Tcut is the cutoff length, which we set 10 fm to cutoff the interaction 
within the length of the cell size. The physical screening length of the Coulomb 
potential by the electron localization is, however, estimated to be much larger 
in the case of the normal nuclear density.'^ Thus our "screening" should be 
considered as a technical approximation to avoid this cell-size dependence and 
to make the numerical calculation feasible. In Sec. 13.21 on the other hand, we 
employ completely uniform electron background and fully include the long-range 
Coulomb interaction of protons in replica cells by the Ewald summation method. 
In both the cases we do not simulate electrons explicitly by QMD. 

§3. Nuclear matter at subsaturation densities by QMD 

3.1. Appearance of inhomogeneous structures in nuclear matter 

Through a number of works, it has been clarified that the pasta phases might 
be the ground state of matter at subsaturation densities for various nuclear inter- 
actions including typical oneS.E3'E2J.!MI,|25}.271,23),281,291,30),31),321,33l,34|,35) 

addition, the pasta phases can occupy a significant mass fraction of neutron star 
crusts (~ 50%^P3 and collapsing supernova cores (> 20%)"^^' if they really exist 
in these objects. 

However, almost all the previous works assume several possible shapes of nuclei 
and what they can actually claim is that the pasta phases can be the energetically 
most favorable state among the selected specific structures. Furthermore, all 
previous studies are based on static frameworks and focus only on the equilibrium 
state, mainly the ground state. Therefore, the fundamental problem whether or 
not the pasta phases are actually formed in young neutron stars in their cooling 
process and supernova cores in the stage of the gravitational collapse has been 
totally unclear. 

To solve this problem, we have studied whether the pasta phases are formed by 
adiabatically changing an external parameter (either decreasing the temperature 
or increasing the density) without any assumption on the nuclear shape.'^''^''22J'E2l.llU.ll2li 
For this purpose, QMD which enables us to simulate the time evolution of the 
nucleon many-body systems with a large number of nucleons is very powerful. 
It is also noted that we are mainly interested in the nuclear structure from the 
mesoscopic to macroscopic scales of > 10 fm, where the exchange effect should 
be less important. Therefore, it is expected that QMD is a reasonable approx- 
imation for studying the pasta phases. Especially, at non-zero temperatures of 
> 0(1) MeV, validity of QMD is ensured because the shell effects are washed out 
by thermal fluctuations above T - 3 MeV.I^S 

We consider a system with neutrons, protons, and electrons in a cubic box 
with periodic boundary conditions. The system is not magnetically polarized, 
i.e., it contains equal numbers of protons (and neutrons) with spin up and spin 
down. Relativistic degenerate electrons which ensure charge neutrality can be 
regarded as a uniform background because electron screening is negligible at rele- 
vant densities around the normal nuclear densitj'^''22J'l2U as we have discussed in 
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the last section. If we assume completely uniform electron distribution, we have 
to take account of the long-range nature of the Coulomb interaction. First, we 
show the results of our QMD calculations with a cutoff distance of the Coulomb 
interaction, in order that the range of the interaction does not exceed the cell 
size.l^ Similar calculation has been done also in a pioneering work by Peilert et 
al. using a different QMD model.® 

Figure [7] is snapshots of symmetric nuclear matter at various densities. Below 
0.6po, the ground state of matter becomes inhomogeneous.'^ We can see struc- 
tures similar to nuclear "pasta". However, these structures are not regular and 
some of them are hard to classify into typical pasta structures. In addition, we 
could not realize the bcc lattice of spherical nuclei at low densities of ^ O.lpo, 



which must be the true ground state in this density region. We consider the 
reasons of these incompleteness and why we could not obtain the pasta phases 
are the use of a cutoff distance for the Coulomb interaction and possibly due to 
too rapid cooling. 

Energy per nucleon of symmetric nuclear matter is plotted by open squares 
in Fig. [6l The contribution of electrons is subtracted and only the nucleon con- 
tribution is included. Compared with uniform distribution (filled squares), the 
energy gets lower at densities below ^ O.Gpo- This is due to the formation of 
inhomogeneous structures. 

Figure [8] shows the dependence of the structure on the proton fraction. If the 
proton fraction x is close to 0.5, neutrons and protons distribute congruently.-^ 
However, with decrease of x, extra neutrons drip out of the nuclei and spread the 
whole space. 

3.2. Pasta formation by decreasing temperature 

The above QMD results confirm that inhomogeneous structures appear at 
subsaturation densities and this phenomenon has a significant effect on the EOS. 
However, an important question whether or not the pasta phases actually ap- 
pear is still unclear. To solve this problem, we have performed another series of 
QMD studie^^''^''^''^ in which we include the long-range contribution of the 
Coulomb interaction and perform careful simulated annealing to achieve thermal 
equilibrium. Here we calculate the Coulomb interaction by the Ewald summation 
method, which enables us to sum up the contributions of long-range interactions 
in a system with periodic boundary conditions efficiently. For nuclear interac- 
tion, we use the QMD Hamiltonian of Ref. 1^0)) with the standard medium EOS 
parameter set and another form of the QMD Hamiltonian of Refs. H7)) .H5 |) . 
The qualitative results are the same for both the models.® 

In Fig.|9l we show the resulting snapshots of the nucleon distributions for x — 
0.3 at r ~ MeV. Note that we obtain the pasta structures without assuming the 
nuclear shapes a prioriW^^^ This is the first result which shows the formation 
of the pasta phases using dynamical framework. In the simulations, we first 
prepare a uniform hot nucleon gas at T 20 MeV for each density. Starting from 



*■* Recently, Horowitz and his collaborators have also studied the structure of nuclear matter at 
subsaturation densities usine 

QMD .1^ 'El) Their model is close to the early version of QMD without 
the Pauh potential, and thus they cannot simulate systems at zero temperature. 

**' Very recently, one of the present authors and his collaborators have found that fee lattice of 
spherical nuclei can be the ground state, by taking the optimum sizes of the cell and nuclei as well 
as the inhomogeneous electron distribution.!^ 





24fm 



p=0.2po 



p=0.1po 




35fm 




58fm 



Fig. 7. Nucleon distributions of symmetric nuclear matter x = 0.5 at T ~ 0. The total number 
of nucleons in the simulation is 1024. The red particles show protons and the white ones show 
neutrons. For densities above O.Spo, matter is uniform. At lower densities, there appear some 
incomplete pasta-like structures: spherical bubbles (0.6po), rod-like nuclei (0.2po), and spherical 
nuclei (O.l/Oo). This figure is taken from Ref. I20p. 
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Z/A=0.293 Z/A=0.391 




;. 8. Nucleon distributions of asymmetric nuclear matter for p = O.lpo at T ~ 0. The dependence 
on the proton fraction x = Z/A is shown. This figure is taken from Ref. I20p . 

this initial condition, we slowly cool it down using frictional relaxation method, 
which is given by QMD equations of motion plus small friction terms [Eq. (|12p ]. 
Throughout this cooling process, we keep the quasi-thermal equilibrium. We take 
the time scale of 0(10'' — 10'*) fm/c to reduce the temperature down to ^ 0.1 MeV 
or less. This result suggests that the pasta phases can be formed in the neutron 
star crusts by cooling. It should be noted, however, that the typical value of the 
proton fraction x in the relevant region of neutron star crusts is < 0.1, which is 
lower than that used in these simulations. However, we have also done similar 
simulations at x = 0.1 and have observed the formation of the pasta phases.'^Sl' 

Another important advantage of QMD is that effects of non-zero temperature 
can be naturally incorporated. Therefore, QMD simulations provide us a clear 
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picture how the pasta structures shown in Fig. [9] are formed by decreasing tem- 
peraturePJ In Figs. [TO] and HH we show snapshots of the nucleon distributions 
at non-zero temperatures for x = 0.3 at p = 0.175po and 0.34/9o, respectively; 
at T ~ MeV, we obtain the pasta phase with rod-hke nuclei in the former 
case, and that with slab-like nuclei in the latter case. Around these densities, 
the phase separation occurs at T ~ 5 MeV and we see that, at T ~ 3 MeV, the 
density inhomogeneity by clustering of nucleons becomes significant [Figs. ITOT c) 
and [TlTc')]. At T ~ 2 MeV, nuclear shapes become recognizable even though 
the surface diffuseness of nuclei and the fluctuation of the nuclear shape are still 
large and there are still many evaporated nucleons among nuclei [Figs. llOf b) and 
ITlT b)]. By further decreasing temperature, these surface diffuseness, fluctuation 
of the nuclear shape, and the number of evaporated nucleons except for dripped 
neutrons become small and, eventually, clear pasta structures can be observed at 
T < 1 MeV [Figs. [TOta) and[TT|;a)]. 

Finally, we summarize our results in a phase diagram at subsaturation densi- 
ties shown in Fig. [T^ In the region below the thick dotted lines, where we can 
identify the nuclear surface, we have obtained the pasta phases with spherical 
nuclei [region (a)], rod-like (cylindrical) nuclei [region (b)], slab-like nuclei [region 
(d)], cylindrical holes [region (f)] and spherical holes [region (g)]. It is noted 
that, in addition to the pasta phases of these simple structures, phases with 
more complicated structures whose both the nuclear matter region and the bub- 
ble region have multiply-connected configurations have been obtained [regions (c) 
and (e)]. Existence of phases with such complicated structures, e.g., gyroid and 
double-diamond phases, have been discussed by several authors using different 
methods ^ ^ ^ ^ ^ 

3.3. Pasta formation by compression of nuclear matter 

In supernova cores, pasta phases are expected to be formed by compression 
of matter in the gravitational collapse [see, e.g., Refs. [5T]),[521) and references 
therein]. This issue is, however, more non-trivial compared to the pasta formation 
by decreasing temperature discussed in Sec. 13.21 because drastic changes of the 
nuclear structure, such as from sphere to rod, must be involved in the present 
case. Therefore, to solve whether or not the pasta phases are formed in collapsing 
supernova cores and to understand its formation process, an ab-initio dynamical 
approach is needed. Using QMD, which is very suitable for this purpose, we have 
solved this problem by demonstrating that a lattice of rod-like nuclei is formed 
from a bcc lattice of spherical nuclei by compression!— 1* ^1. Our results establish 
that the pasta phases can be formed in collapsing supernova cores. 

A generally accepted scenario of the formation of the pasta phases in supernova 
cores is that, when the density exceeds some critical value, an instability of nuclear 
fission sets in and, consequently, all the nuclei elongate in the same direction and 
"eventually join up to form string-like structures" ~^ When the Coulomb energy 
between protons in a nucleus is sufficiently larger than the surface energy of 
the nucleus, the reduction of the Coulomb energy due to the fission exceeds the 
energy cost of the surface tension by an increase of the surface area, and the 
fission barrier vanishes (i.e., the onset of the fission instability). This is described 



*' We have shown that a layered lattice of slab-like nuclei and triangular lattice of cylindrical 
bubbles are also formed by further compressiorP^ 
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Fig. 9. Nucleon distributions of tiie pasta piiases for x = 0.3 at T ~ MeV. The total number 
of nucleons in the simulation is 2048 (614 protons and 1434 neutrons). The red particles show 
protons and the green ones show neutrons. Each panel shows the pasta phase with (a) spherical 
nuclei (O.lpo), (b-1) rod-like nuclei (O.lSpo); side view, (b-2) the same; top view, (c) slab-like 
nuclei (0.35po), (d) rod-like bubbles (0.5po), and (e) spherical bubbles (0.55po)- This figure is 
adapted from Ref . I38p . 
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(a)r=1MeV (b)r=2MeV (c)r=3MeV 




Fig. 10. Nucleon distributions at T = 1, 2, and 3 MeV for x = 0.3 and p — 0.175po, where the 
phase with rod-hke nuclei is obtained at zero temperature. The total number of nucleons in this 
simulation is 16384 (4915 protons and 11469 neutrons). The upper figures show top views along 
the axis of the rod-like nuclei at T = 0, and the lower ones show side views. The red particles 
represent protons and green ones represent neutrons. This figure is taken from Ref. I39|) . 



(a)r=1MeV (b)r=2MeV (c)r=3MeV 




Fig. 11. The same as Fig. [TO] for x — 0.3 and p = 0.34po, where the phase with slab-like nuclei is 
obtained at zero temperature. The line of sight of these figures is in the direction parallel to 
the plane of the slab-like nuclei at T = 0. This figure is taken from Ref. I39p . 

by the celebrated Bohr- Wheeler condition: 

-Bcoul > SE'surf, (21) 

where E^^^^^ and -Esurf are the self-Coulomb and surface energies of the nucleus. 
Using an equilibrium condition between the Coulomb and surface energies eval- 
uated within the WS approximation, one can show that Eq. (I^TI) reads u > 1/8, 
where u is the volume fraction occupied by nuclei.'^ 

However, we should note that Bohr- Wheeler condition (PT|) is derived for an 
isolated nucleus in vacuum. In the real situation in supernova cores, there are 
background electrons and the condition for the fission instability should be mod- 
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Fig. 12. Phase diagram of matter at a; = 0.3 plotted in the p-T plane. The dashed lines correspond 
to phase separation lines. The thick dotted lines show the boundary above which nuclear 
surface cannot be identified. The dashed-dotted lines show boundaries between different phases. 
Abbreviations SP, C, S, CH, and SH mean phases with spherical nuclei, cylindrical (rod-like) 
nuclei, slab-like nuclei, cylindrical holes, and spherical holes, respectively. The parentheses (A,B) 
show an intermediate phase between A and B phases. Simulations have been carried out at the 
points denoted by circles. This figure is adapted from Ref. |39j). 



ified from the original Bohr- Wheeler condition ([2T|). Indeed, it has been shown 
that the fission instability is suppressed by the background electrons which reduce 
the local net charge density inside nuclei.'^''^ This result poses a doubt about 
the formation scenario based on the fission instability. 

In Fig. I13[ we show the snapshots of our simulation, which capture the forma- 
tion of the pasta phase in adiabatic compression. Here we use the QMD Hamilto- 
nian of Ref. I20p with the standard medium EOS parameter set as in the previous 
section. The proton fraction x and the total number of particles N are x ~ 0.39 
and N — 3328 (with 1312 protons and 2016 neutrons) in this simulation. Starting 
from an initial condition at p = O.lSpo and T = 0.25 MeV [Fig.fT^Ta)]. we increase 
the density by changing the box size Lceii slowly at a rate < O(10~^) po/{fm/c), 
which yields the time scale of > 10^ fm/c to reach the typical density region of 
the phase with rod-Uke nuclei r^l At p ~ 0.243po [Fig. [T^T c)]. the first pair of two 
nearest-neighbor nuclei start to touch and fuse (dotted circle), and then form an 



*' While this time scale is, of course, much smaller than the actual time scale of the collapse, it is 
much larger than that of the change of the nuclear shape (e.g., ~ 1000 fm/c for the nuclear fission) 
and thus the dynamics observed in the simulation should be governed by the intrinsic physical 
properties of the system, not by the density change applied externally. 
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elongated nucleus. Then, multiple pairs of nuclei fuse and become such elongated 
nuclei in a way that they are aligned in a zigzag configuration [Fig.fTSfd)]. These 
elongated nuclei further stick together [see Figs. [T^e) and (f)], and all the nu- 
clei fuse to form rod-like nuclei as shown in Fig. fTHTg). Finally, we obtain an 
almost perfect triangular lattice of rod-like nuclei after relaxation [Figs. IT^ h-l) 
and (h-2)]. 

Note that, until the nuclei touch and fuse, they keep the spherical shape [see 
Fig. ITST c)]. This shows that the pasta phase is formed without undergoing the 
fission instability. It is also remarkable that, in the middle of the transition 
process, the elongated nuclei made of a pair of spherical nuclei take a zigzag 
configuration and then they further connect to form wavy rod-like nuclei. This 
process is very different from the above-mentioned conjectured scenario based on 
the fission instability. 

3.4. Expansion of nuclear matter 

In this section, we present another example of a MD simulation of matter, 
i.e., expanding nuclear matter .l^'l^ Multifragmentation is one of the topics of a 
long-standing interest in heavy-ion collision physics. Not only does its mechanism 
itself attract our interest but also it is a good test bed for the EOS. In particular, 
fragment mass spectra have been discussed in many works. The fragmentation 
mechanism is different between the participant region where colliding two nuclei 
overlap, and the spectator region which surrounds the participant region. In the 
participant region, the fragmentation is affected by a radial flow, and thus this 
region is not in thermal equilibrium. Empirically, a large radial flow yields an 
exponential mass spectra.'^''^ On the other hand, for the spectator region, 
statistical models are effective since the collective dynamics is not important. 
Fisher's droplet model, which is one of the most famous statistical models, predict 
a power law of fragment mass spectra. 

There are two major scenarios for fragment formation. One is as follows: by 
the collision of two nuclei, a hot and dense zone is made. It expands due to 
the high temperature and high densities. At this point, the system is considered 
as a gas. As the system becomes dilute, the temperature decreases and finally 
crosses the liquid-gas transition point. Then the system experiences the spinodal 
instability and fragmentation occurs. The other scenario is that the expanding 
region behaves like a solid. Then the fragmentation occurs by the formation of 
cracks. 

Many studies have been performed so far for coflision analysiiP^ -l^'l^'lsSJ .113.111.1131. 113 
which include both dynamics and statistics. However, the dynamics are compli- 
cated due to the finite size effect. In order to simplify the problem and to get 
more direct insights into the fragmentation mechanism, we perform QMD simula- 
tions of infinite matter by imposing special periodic boundary conditions with an 
isotropic expansion.!^ Prior to our study, instability of nuclear matter against 
multifragmentation has been studiecp3 and also a simulation study of the infinite 
system based on a dynamical model has been tried by other authors. These 
studies, however, do not take into account the expansion of the system. 

3.4.1. Model 

In order to simulate expanding nuclear matter, we employ a special periodic 
boundary condition. A method to simulate two-dimensional (2D) expanding pe- 
riodic systems has been proposed in condensed matter physics. fi^J .EH We extend 
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Fig. 13. Snapshots of the formation process of the pasta phase with rod-hke nuclei from a bcc 
lattice of spherical nuclei by compression of matter. The red particles show protons and the 
green ones show neutrons. In panels (a)-(g) and (h-1), nucleons in a limited region [surrounded 
by the dotted lines in panel (h-2)] are shown for visibility. The vertices of the dashed lines in 
panels (a) and (d) show the equilibrium positions of nuclei in the bcc lattice and their positions 
in the direction of the line of sight are indicated by the size of the circles: vertices with a large 
circle, with a small circle, and those without a circle are in the first, second, and third lattice 
plane, respectively. The solid lines in panel (d) represent the direction of the two elongated 
nuclei: they take zigzag configuration. The box sizes are rescaled to be equal in the figures. 
This figure is taken from Ref. I42[>. 



this method for 3D systems and apply it to nuclear fragmentation.SS-BS) 

Hamiltonian used here consists of the following effective interaction terms. 



H = J2iif+ "^""=1 + ^-'-f + ^P-"ii' (22) 
1=1 
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where (p^) = X^j^^^^i) Pij defined in Eqs. (O and (O. Although this Hamiho- 
nian is almost the same as that used in the previous sections (introduced in Sec. 
12. 2p . a surface term Vsmf is newly added and parameters are readjusted as shown 
in Table HH By this improvement of the effective interaction, experimental data 
of the radii and the binding energies of nuclei are reproduced better. Coulomb 
potential and a two-body collision term are not included for simplicity. 



Table II. Effective interaction parameter set (_fsr=280 MeV). 



a (MeV) 


-121.9 


13 (MeV) 


197.3 


7 


4/3 


Cso (MeV) 


25.0 


(MeV) 


-258.5 


Cex(2) (MeV) 


375.6 


Ml (MeV) 


2.35 


/i2 (MeV) 


0.4 


VsF (MeV) 


20.68 


Cp (MeV) 


115.0 


PO (MeV) 


120.0 


go (fm) 


2.5 


A (fm) 


1.40 







Numerical calculations are done in the primitive cubic cell and 26 replica cells 
surrounding the primitive cell. With the standard periodic boundary condition, 
we first prepare a set of equilibrated initial conditions of symmetric nuclear mat- 
ter at Po for various values of the temperature. In this model, we set the nuclear 
saturation density po = 0.16 fm~^. According to the Metropolis sampling proce- 
dure, 1000 samples are created for each temperature. Next, we give each nucleon 
an additional collective momentum proportional to its position vector. The same 
motion is given to each cell boundary in order that the whole system undergoes 
a homogeneous expansion. The collective momentum PcoU of a particle i in the 
replica cell is also proportional to the position vector R^ -I- Lccii of the particle [see 
Fig. [T3] . The speed of expansion is characterized by a radial flow velocity param- 
eter ft,, which is analogous to the Hubble constant in cosmology. The collective 
momentum added to a particle at a position R is as follows: 

Peoll(R) = h^^Pv, (26) 
Po 

where Pp is the Fermi momentum at saturation density po- Each sample with a 
given temperature will be expanded with radial flow velocity parameter h. Note 
that when a nucleon goes out of a boundary of the primitive cell, its image particle 
comes in from the opposite boundary of the cell. What is different from the normal 
periodic boundary conditions is that the momentum of the image particle is also 
modified by the collective momentum Pcoii(Lceii) proportional to the cell size. 
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Fig. 14. Schematic explanation of the periodic boundary condition with uniform expansions. The 
left box is the primitive cell and the right one is a replica cell located at Lccii. The cell velocity 
is Vceii = Lcell. This figure is taken from Ref. I49|> . 

In our setup, this collective motion does not stop. Moreover, it does not 
change at all. This is because the collective force per unit cell which intends 
to change the collective momentum is finite but the moment of inertia per unit 
cell diverges for the infinite system. In other words, while the potential energy 
per cell is finite, the kinetic energy per cell is infinite due to the contributions of 
infinitely distant cells with diverging velocities. We cannot avoid this constant 
collective velocity when we consider expansion or contraction of infinite systems 
with periodic boundary conditions!*^ 




Fig. 15. Snapshots of nucleon distribution in the primitive cell in the initial state at po (left panel) 
and the final state at 0.05po (right panel). This figure is adapted from Ref. I49[>. 

We calculate the time evolution until the average density reaches O.OSpo, below 
which we can identify fragments. Figure [15] shows an example of our simulation; 
the initial state at po and T = 30 MeV (left panel) and the final state at O.OSpo 
(right panel). Each fragment is determined by the clustering algorithm: When 
the distance between a particle and a fragment is smaller than 3 fm, we identify 
that this particle belongs to the fragment. We perform 1000 events for each h 
and T to get fragment mass spectra. 



*•* In the simulation of adiabatic compression of nuclear matter (Sec. 13.3]) . we have rescaled the 
particle positions at the compression rate instead of using the initial collective momentum explained 
in this section. This is valid because the compression of the system in that case is very slow compared 
to the motion of nucleons. 
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3.4.2. Results 

Before discussing how fragment mass spectrum depends on h and T, we study 
how the pressure of the system is related to the dynamical expansion. The in- 
stantaneous pressure is calculated on the basis of the virial theorem as follows!^ 

O "pthermal i 

P - i^PE + 6^^5:E(I^' - I^.) • (27) 

Z I j 

where = is the velocity of particle i and F.^ is the force acting between 
particles i and j. Note that the momentum includes two components: One is the 
momentum of the thermal motion, and the other is that of the collective motion 
characterized by h. To calculate the pressure, we subtract the momentum of the 
collective component of Eq. from P„ i.e., pthormai ^ _ pcoU 




0.5 1.0 
P/Po 



Fig. 16. Upper panels: Density dependence of the pressure. From the left, the initial temperature 
Tlnit = 30, 15, and 5 MeV. Dots show isotherms which are calculated at each fixed density by 
the Metropolis sampling method. Lines show adiabatic cases that nuclear matter expands from 
saturation density po to O.OSpo. Lower panels: Density dependence of the effective temperature. 
This figure is taken from Ref. 149 [I. 

Figure \W\ shows the relation between density and pressure. Bold lines are 
isothermal (T = 5, 15, and 30 MeV) cases calculated by the Metropolis sam- 
pling method at each density. Thin lines are the trajectories along which the 
system expands. The instantaneous pressure of the expanded matter is smaller 
than that of the equilibrium state. One reason of this reduction is that the ex- 
pansion of matter proceeds adiabatically, i.e., without an exchange of heat and 
with a positive work to the environment. Therefore, the system is cooled down 
as it expands Q Another reason is that the expansion is a dynamical and non- 



In the case of low initial temperatures, there appears a region of negative effective tem- 
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Fig. 17. The fragment mass spectra at p = 0.05po- The fragment muhiphcity M is counted in 
the primitive cell. The expansion begins at saturation density po with the initial configuration 
created by Metropolis sampling method. Upper figure shows the results for initial temperature 
5 MeV and lower one for 30 MeV. The radial fiow velocities corresponding to h = 0.10, 0.50, 
1.00, and 2.00 are shown. Both the right and the left show the same data with different scale 
of abscissa. This figure is taken from Ref. I48p . 



equilibrium process. If the system were static, the region of mechanical instability 
with a negative gradient of the pressure against the density would be avoided by 
a formation of clusters0 For the expanding case, on the other hand, the system 
cannot reach the stable state with a higher pressure. This is the same tendency 
observed experimentally in Ref. 1721) . 

The fragment mass distributions are very sensitive to the radial flow velocity. 
Figure [T7] shows the fragment mass distributions obtained at p = O.OSpo after 
the dynamical expansion. Both the right and the left panels show the same data 
but with different scales of abscissa. As the radial flow velocity increases, the 
fragment mass distribution changes from the power law, i.e., a straight line in the 

peratures (see the lower right panel of Fig. llOp . This means that some particles have velocity 
and momentum of opposite directions, which is caused by the momentum-dependent terms of the 
potential. 

*' One may notice that the pressure in the upper right panel of Fig. [16] still has a negative 
gradient. However, with inclusion of electron contribution, the total pressure shows monotonically 
increasing dependence on the density. Similar situation is reported in Ref. I7ip . 
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double logarithmic scale, to the exponential decay, i.e., a straight line in the linear- 
logarithmic scale. It is also noted that the fragment mass distribution decreases 
more rapidly with increasing h. This feature is consistent with experimental data 
of heavy-ion collisions at small impact parameters where a large radial flow is 
observed .'SS 

In the fast expansion limit, the effect of the interaction between particles is 
expected to be small and the fragmentation is determined solely by the initial 
position of particles. The fragment mass distribution in this limit shows an ex- 
ponential forrcPSJ and the change of the slope according to the expansion velocity 
is also consistent with the consideration given in Ref . f62j) . 

In the case of slow expansion {h = 0.1), the fragment mass distributions 
for Tinit — 5 and 30 MeV agree very well to the Fisher's power law. One may 
conclude that these fragment mass distributions indicate an evidence of the liquid- 
gas phase transition of nuclear matter. However, these results do not always 
support this scenario. The reason is as follows: One important premise of the 
Fisher's model is that the temperature is above the critical temperature, which 
is about 8 MeV in the present model.!^ On the other hand, in our case with the 
initial temperature of 5 MeV, where the mass distribution shows a power law, the 
freeze-out temperature should be lower than 5 MeV. Therefore, the fragment mass 
distribution with the power low is not necessarily accompanied by the liquid-gas 
phase transition. Instead, a breakup of a solicP^'l^ by cracking is more plausible 
in the case of lower temperatures. 

§4. Summary 

In this article, we have reviewed the molecular dynamics approaches to nu- 
cleon many-body systems, focusing on the quantum molecular dynamics (QMD) 
model. This method was originally developed for the study of heavy-ion collisions 
to describe multifragmentation which the time-dependent Hartree-Fock (TDHF) 
theory failed to explain. With the help of the Pauli potential to take account of 
the Pauli principle, QMD has been applied to dense nuclear matter in addition 
to heavy-ion collisions. 

A great advantage of MD approaches is that we can study the dynamical 
process of nucleon many-body systems without any assumptions on the nuclear 
structure. QMD, in particular, is very suitable to study inhomogeneous nuclear 
matter in the mesoscopic to macroscopic scales because QMD makes it possible 
to simulate large systems with many nucleons. To show the power of this method, 
we have presented some of our works using QMD. 

First, we have explained our QMD model in Sec. [2j We have seen that this 
model is designed to give the correct saturation properties and reasonable equa- 
tions of state (EOS) of nuclear matter, and give a good agreement of the binding 
energies of light nuclei with ^ ^ 8 including alpha particle and also heavy nuclei 
with A > 40. These points are important for the reliability of the predictions by 
this model, which we have presented in the remaining part of this article. 

From Sec. l3.1l to l3.31 we have shown a series of our studies about the structure 
of nuclear matter at subsaturation densities. It has been predicted that nuclear 
"pasta" phases, states of matter with nuclei of rod-like and slab- like shapes, can 
be the ground state of matter in this density region. Using QMD, we have shown 
that the pasta phases can actually be formed by cooling hot uniform nuclear 
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matter and by compressing a bee lattice of spherical nuclei. These results strongly 
suggest that the pasta phases are formed in the cooling process of hot neutron 
star crusts and by the compression of matter in the collapse of supernova cores. 

These results have important implications on the mechanical strength of the 
neutron star crust, the cooling process of hot protoneutron stars, the mechanism 
of glitches, etc. Making an EOS table for core collapse simulations, taking into 
account the existence of the pasta phases including their effects on the neutrino 
opacity, is an important direction. 

In Sec. 13.41 we have discussed fragment formation in expanding nuclear matter. 
We have developed a method to describe isotropically expanding matter using 
the periodic boundary conditions. Using this method together with our QMD 
model, we have calculated the fragment mass distribution and have found that it 
shows an exponential decay for rapid expansion and a power-law decay for slow 
expansion. Our analysis suggests that multifragmentation at lower temperatures 
occurs not by the liquid-gas phase transition but by some other mechanisms, e.g., 
the formation of cracks in the solid-like expanding region. 

Molecular dynamics simulations have been playing important roles in nuclear 
physics: both in the study of the nuclear structure and reaction. They can 
successfully describe the structure of nuclei and statistical properties of heavy-ion 
collisions taking account of many-body correlations and fluctuations. In addition 
to nucleon many-body systems, MD simulations are used also in QCD studies: 
UrQMD is now commonly used to analyze heavy-ion collisions with a quark-gluon 
plasma,'^! and some dynamical properties of quark matter have been studied by 
a MD approach.ESJ''^ One of the most important and challenging directions is 
to incorporate the wave nature of the quantum mechanics in the MD approach, 
which is based on the particle picture. A new framework beyond QMD and 
FMD/AMD in this direction is highly awaited.'^ By such a breakthrough, MD 
should be a promising approach also for studying the dynamics of fission and 
fusion, which is a long-standing mportant problem in nuclear physics. 
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